Load libraries

library(readxl)
library(tidyverse)
library(jyluMisc)
library(DrugScreenExplorer)
library(ComplexHeatmap)
library(ggrepel)
library(parallel)
library(gridExtra)
library(knitr)

Define variables

opt <- list()
opt$drugscreen <- "data/submission/drugScreens_pseudo.RDS"
opt$druganno <- "misc/drugList_suppl.xlsx"
opt$plot <- "plots/SFig2/"

## Set theme
lgd <-  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5), 
        text = element_text(size = 15),
        panel.background = element_rect(fill = "transparent",colour = NA),
        plot.background = element_rect(fill = "transparent",colour = NA), 
        legend.background = element_rect(fill='transparent',colour = NA), 
        strip.background = element_blank(), 
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

Load drug screen data and annotations

drugList <- readRDS(opt$drugscreen)
drugAnno <- read_excel(opt$druganno) %>%
  dplyr::select(-Supplier, -Screen) %>% unique()

Analysis

Dose-response verdinexor

embl2016 <- drugList[["ScreenB"]]

drugSub <- "Verdinexor"

concInter <- embl2016 %>%
  dplyr::filter(name %in% c(drugSub), 
         diagnosis == "T-PLL") %>%
  dplyr::select(conc) %>% unique() %>% unlist() %>% 
  as.vector()
concInter[1:2] <- round(concInter[1:2], 1)
concInter[3:6] <- round(concInter[3:6], 2)
concInter[7:9] <- round(concInter[7:9], 3)

p1 <- embl2016 %>% 
  dplyr::filter(name %in% c(drugSub), 
         diagnosis == "T-PLL") %>%
  ggplot(aes(x = conc, y = viab, group = sampleID)) +
  geom_line(size = 0.25) + 
  geom_point(size = 3.75, shape = 21, fill = "#64B5F6") +
  geom_hline(yintercept = 1, linetype = "dashed") +
  scale_y_continuous(limits = c(0, 1.2)) +
  scale_x_log10(labels = concInter, 
                breaks = concInter) +
  xlab("Concentration (µM)") + ylab("Viability") + 
  theme_classic() + ggtitle(paste("Screen B", drugSub, sep = " - ")) +
  theme(text = element_text(size = 17.5), 
        legend.key = element_blank(),
        plot.title = element_text(hjust=0.5),
        # axis.text.x = element_text(hjust = 1, angle = 45),
        axis.text = element_text(size = 17.5),
        axis.line = element_line(linewidth = 0.75),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "transparent",colour = NA),
        plot.background = element_rect(fill = "transparent",colour = NA))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
p1

ggsave(plot = p1, filename = paste0(opt$plot, drugSub, ".png"), 
       height = 5.25, width = 7)

Dose-response ruxolitinib

drugSub <- "Ruxolitinib"

concInter <- embl2016 %>%
  dplyr::filter(name %in% c(drugSub), 
         diagnosis == "T-PLL") %>%
  dplyr::select(conc) %>% unique() %>% unlist() %>% 
  as.vector()
concInter[1:2] <- round(concInter[1:2], 1)
concInter[3:6] <- round(concInter[3:6], 2)
concInter[7:9] <- round(concInter[7:9], 3)

p1 <- embl2016 %>% 
  dplyr::filter(name %in% c(drugSub), 
         diagnosis == "T-PLL") %>%
  ggplot(aes(x = conc, y = viab, group = sampleID)) +
  geom_line(size = 0.25) + 
  geom_point(size = 3.75, shape = 21, fill = "#64B5F6") +
  geom_hline(yintercept = 1, linetype = "dashed") +
  scale_y_continuous(limits = c(0, 1.2)) +
  scale_x_log10(labels = concInter, 
                breaks = concInter) +
  xlab("Concentration (µM)") + ylab("Viability") + 
  theme_classic() + ggtitle(paste("Screen B", drugSub, sep = " - ")) +
  theme(text = element_text(size = 17.5), 
        legend.key = element_blank(),
        plot.title = element_text(hjust=0.5),
        # axis.text.x = element_text(hjust = 1, angle = 45),
        axis.text = element_text(size = 17.5),
        axis.line = element_line(linewidth = 0.75),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "transparent",colour = NA),
        plot.background = element_rect(fill = "transparent",colour = NA))
p1

ggsave(plot = p1, filename = paste0(opt$plot, drugSub, ".png"), 
       height = 5.25, width = 7)

Dose-response duvelisb

drugSub <- "Duvelisib"

concInter <- embl2016 %>%
  dplyr::filter(name %in% c(drugSub), 
         diagnosis == "T-PLL") %>%
  dplyr::select(conc) %>% unique() %>% unlist() %>% 
  as.vector()
concInter[1:2] <- round(concInter[1:2], 1)
concInter[3:6] <- round(concInter[3:6], 2)
concInter[7:9] <- round(concInter[7:9], 3)

p1 <- embl2016 %>% 
  dplyr::filter(name %in% c(drugSub), 
         diagnosis == "T-PLL") %>%
  ggplot(aes(x = conc, y = viab, group = sampleID)) +
  geom_line(size = 0.25) + 
  geom_point(size = 3.75, shape = 21, fill = "#64B5F6") +
  geom_hline(yintercept = 1, linetype = "dashed") +
  scale_y_continuous(limits = c(0, 1.2)) +
  scale_x_log10(labels = concInter, 
                breaks = concInter) +
  xlab("Concentration (µM)") + ylab("Viability") + 
  theme_classic() + ggtitle(paste("Screen B", drugSub, sep = " - ")) +
  theme(text = element_text(size = 17.5), 
        legend.key = element_blank(),
        plot.title = element_text(hjust=0.5),
        # axis.text.x = element_text(hjust = 1, angle = 45),
        axis.text = element_text(size = 17.5),
        axis.line = element_line(linewidth = 0.75),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "transparent",colour = NA),
        plot.background = element_rect(fill = "transparent",colour = NA))
p1

ggsave(plot = p1, filename = paste0(opt$plot, drugSub, ".png"), 
       height = 5.25, width = 7)

Show reproducibility

diagList <- c("T-PLL")

screenNames <- c("ScreenB", "ScreenC", "ScreenD", "ScreenA", "ScreenE")

pList <- lapply(diagList, function(x) {
  pList <- lapply(screenNames, function(y) {
    screenData <- drugList[[y]]
    
    screenData <- screenData %>% #filter(diagnosis == x) %>%
      group_by(diagnosis, patientID, name) %>%
      summarise(viab.auc = mean(normVal, na.rm = TRUE)) %>%
      ungroup()
    
    screenNamesOther <- screenNames[which(screenNames != y)]
    
    pList <- lapply(screenNamesOther, function(z) {
      screenDataOther <- drugList[[z]]
    
      screenDataOther <- screenDataOther %>% #filter(diagnosis == x) %>%
        group_by(diagnosis, patientID, name) %>%
        summarise(viab.auc.other = mean(normVal, na.rm = TRUE)) %>%
        ungroup()
      
      ## Intersect the patients
      patOver <- intersect(screenData$patientID, screenDataOther$patientID) %>% unique()
      message("Found ", length(patOver), " overlapping patients")
      
      ## Intersect the drugs
      drugOver <- intersect(screenData$name, screenDataOther$name) %>% unique()
      message("Found ", length(drugOver), " overlapping drugs")
      
      ## Plot overlap
      p <- screenData %>%
        filter(patientID %in% patOver, name %in% drugOver) %>%
        left_join(screenDataOther, by = c("name", "patientID", "diagnosis")) %>%
        ggplot(aes(x = viab.auc, y = viab.auc.other, fill = diagnosis)) +
        geom_point(size = 2, shape = 21) +
        geom_smooth(method = "lm", colour = "black", fill = "grey80") +
        xlab(paste0(gsub(x = y, pattern = "Screen", replacement = "Screen "))) +
        ylab(paste0(gsub(x = z, pattern = "Screen", replacement = "Screen "))) +
        xlim(0, 1.6) + ylim(0, 1.6) +
        #ggtitle(paste0("Drug-Drug corr. ", y, " vs ", z)) +
        lgd +
        theme(legend.position = "none", 
              panel.border = element_rect(linewidth = 1, fill = "transparent")) +
        scale_fill_manual(values = c("MCL" = "#F44336", "T-PLL" = "#64B5F6", "CLL" = "#FFD45F"))
      
      ggsave(plot = p, filename = paste0(opt$plot, "/screen_corr/", x, "_", y, "_", z, ".png"), 
             height = 7, width = 7.5, units = "cm")
      p
    })
    #grid.arrange(grobs = pList, ncol = 1)
    names(pList) <- screenNamesOther
    pList
  }) 
  names(pList) <- screenNames
  pList
})
names(pList) <- diagList

pList
## $`T-PLL`
## $`T-PLL`$ScreenB
## $`T-PLL`$ScreenB$ScreenC

## 
## $`T-PLL`$ScreenB$ScreenD

## 
## $`T-PLL`$ScreenB$ScreenA

## 
## $`T-PLL`$ScreenB$ScreenE

## 
## 
## $`T-PLL`$ScreenC
## $`T-PLL`$ScreenC$ScreenB

## 
## $`T-PLL`$ScreenC$ScreenD

## 
## $`T-PLL`$ScreenC$ScreenA

## 
## $`T-PLL`$ScreenC$ScreenE

## 
## 
## $`T-PLL`$ScreenD
## $`T-PLL`$ScreenD$ScreenB

## 
## $`T-PLL`$ScreenD$ScreenC

## 
## $`T-PLL`$ScreenD$ScreenA

## 
## $`T-PLL`$ScreenD$ScreenE

## 
## 
## $`T-PLL`$ScreenA
## $`T-PLL`$ScreenA$ScreenB

## 
## $`T-PLL`$ScreenA$ScreenC

## 
## $`T-PLL`$ScreenA$ScreenD

## 
## $`T-PLL`$ScreenA$ScreenE

## 
## 
## $`T-PLL`$ScreenE
## $`T-PLL`$ScreenE$ScreenB

## 
## $`T-PLL`$ScreenE$ScreenC

## 
## $`T-PLL`$ScreenE$ScreenD

## 
## $`T-PLL`$ScreenE$ScreenA

screenNames
## [1] "ScreenB" "ScreenC" "ScreenD" "ScreenA" "ScreenE"
## Give R2 and p.adj
cor.df.all <- lapply(diagList, function(x) {
  lapply(screenNames, function(y) {
    screenData <- drugList[[y]]
    
    screenData <- screenData %>% #filter(diagnosis == x) %>%
      group_by(diagnosis, patientID, name) %>%
      summarise(viab.auc = mean(normVal, na.rm = TRUE)) %>%
      ungroup()
    
    screenNamesOther <- screenNames[which(screenNames != y)]
    
    lapply(screenNames, function(z) {
      screenDataOther <- drugList[[z]]
    
      screenDataOther <- screenDataOther %>% #filter(diagnosis == x) %>%
        group_by(diagnosis, patientID, name) %>%
        summarise(viab.auc.other = mean(normVal, na.rm = TRUE)) %>%
        ungroup()
      
      ## Intersect the patients
      patOver <- intersect(screenData$patientID, screenDataOther$patientID) %>% unique()
      message("Found ", length(patOver), " overlapping patients")
      
      ## Intersect the drugs
      drugOver <- intersect(screenData$name, screenDataOther$name) %>% unique()
      message("Found ", length(drugOver), " overlapping drugs")
      
      ## Plot overlap
      cor.df <- screenData %>%
        filter(patientID %in% patOver, name %in% drugOver) %>%
        left_join(screenDataOther, by = c("name", "patientID", "diagnosis")) %>%
        #group_by(name) %>%
        nest() %>%
        mutate(m = map(data, ~cor.test(~viab.auc+viab.auc.other,.))) %>%
        mutate(res = map(m, broom::tidy)) %>%
        unnest(res) %>% ungroup() %>%
        select(estimate, p.value ) %>%
        arrange(p.value)
      cor.df$screen1 <- y
      cor.df$screen2 <- z
      cor.df
    }) %>% bind_rows()
  }) %>% bind_rows()
}) %>% bind_rows() %>% mutate(p.adj = p.adjust(p.value, method = "BH"))

kable(cor.df.all %>% arrange(p.adj))
estimate p.value screen1 screen2 p.adj
1.0000000 0 ScreenB ScreenB 0
0.7144431 0 ScreenB ScreenC 0
0.6434237 0 ScreenB ScreenD 0
0.7144431 0 ScreenC ScreenB 0
1.0000000 0 ScreenC ScreenC 0
0.7536054 0 ScreenC ScreenD 0
0.6434237 0 ScreenD ScreenB 0
0.7536054 0 ScreenD ScreenC 0
1.0000000 0 ScreenD ScreenD 0
1.0000000 0 ScreenA ScreenA 0
1.0000000 0 ScreenE ScreenE 0
0.5903928 0 ScreenC ScreenA 0
0.5903928 0 ScreenA ScreenC 0
0.6714243 0 ScreenB ScreenA 0
0.6714243 0 ScreenA ScreenB 0
0.8022378 0 ScreenB ScreenE 0
0.8022378 0 ScreenE ScreenB 0
0.8022696 0 ScreenD ScreenE 0
0.8022696 0 ScreenE ScreenD 0
0.6179082 0 ScreenD ScreenA 0
0.6179082 0 ScreenA ScreenD 0
0.8086957 0 ScreenC ScreenE 0
0.8086957 0 ScreenE ScreenC 0
0.6758175 0 ScreenA ScreenE 0
0.6758175 0 ScreenE ScreenA 0
#View(cor.df.all)
max(cor.df.all$p.adj)
## [1] 1.133701e-11
range(cor.df.all[cor.df.all$estimate < 0.99, ]$estimate)
## [1] 0.5903928 0.8086957
## Also per drug
cor.df <- lapply(diagList, function(x) {
  lapply(screenNames, function(y) {
    screenData <- drugList[[y]]
    
    screenData <- screenData %>% #filter(diagnosis == x) %>%
      group_by(diagnosis, patientID, name) %>%
      summarise(viab.auc = mean(normVal, na.rm = TRUE)) %>%
      ungroup()
    
    screenNamesOther <- screenNames[which(screenNames != y)]
    
    lapply(screenNamesOther, function(z) {
      screenDataOther <- drugList[[z]]
    
      screenDataOther <- screenDataOther %>% #filter(diagnosis == x) %>%
        group_by(diagnosis, patientID, name) %>%
        summarise(viab.auc.other = mean(normVal, na.rm = TRUE)) %>%
        ungroup()
      
      ## Intersect the patients
      patOver <- intersect(screenData$patientID, screenDataOther$patientID) %>% unique()
      message("Found ", length(patOver), " overlapping patients")
      
      ## Intersect the drugs
      drugOver <- intersect(screenData$name, screenDataOther$name) %>% unique()
      message("Found ", length(drugOver), " overlapping drugs")
      
      ## Plot overlap
      cor.df <- screenData %>%
        filter(patientID %in% patOver, name %in% drugOver) %>%
        left_join(screenDataOther, by = c("name", "patientID", "diagnosis")) %>%
        group_by(name) %>%
        nest() %>%
        mutate(m = map(data, ~cor.test(~viab.auc+viab.auc.other,.))) %>%
        mutate(res = map(m, broom::tidy)) %>%
        unnest(res) %>% ungroup() %>%
        select(name, estimate, p.value ) %>%
        arrange(p.value) 
      cor.df$screen1 <- y
      cor.df$screen2 <- z
      
      cor.df
    }) %>% bind_rows()
  }) %>% bind_rows()
}) %>% bind_rows() %>% mutate(p.adj = p.adjust(p.value, method = "BH"))
cor.df %>% filter(name %in% c("Birinapant", "Selinexor"))
## # A tibble: 8 × 6
##   name       estimate   p.value screen1 screen2    p.adj
##   <chr>         <dbl>     <dbl> <chr>   <chr>      <dbl>
## 1 Selinexor     0.338 0.000174  ScreenB ScreenD 0.000671
## 2 Birinapant    0.755 0.000120  ScreenB ScreenE 0.000514
## 3 Selinexor     0.306 0.309     ScreenB ScreenE 0.399   
## 4 Selinexor     0.338 0.000174  ScreenD ScreenB 0.000671
## 5 Selinexor     0.912 0.0000355 ScreenD ScreenE 0.000177
## 6 Birinapant    0.755 0.000120  ScreenE ScreenB 0.000514
## 7 Selinexor     0.306 0.309     ScreenE ScreenB 0.399   
## 8 Selinexor     0.912 0.0000355 ScreenE ScreenD 0.000177

Output session info

## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] parallel  grid      stats     graphics  grDevices datasets  utils    
## [8] methods   base     
## 
## other attached packages:
##  [1] knitr_1.48               gridExtra_2.3            ggrepel_0.9.5           
##  [4] ComplexHeatmap_2.18.0    DrugScreenExplorer_0.1.0 jyluMisc_0.1.5          
##  [7] lubridate_1.9.3          forcats_1.0.0            stringr_1.5.1           
## [10] dplyr_1.1.4              purrr_1.0.2              readr_2.1.5             
## [13] tidyr_1.3.1              tibble_3.2.1             ggplot2_3.5.1           
## [16] tidyverse_2.0.0          readxl_1.4.3            
## 
## loaded via a namespace (and not attached):
##   [1] splines_4.3.2               later_1.3.2                
##   [3] bitops_1.0-7                cellranger_1.1.0           
##   [5] lifecycle_1.0.4             Rdpack_2.6                 
##   [7] rstatix_0.7.2               doParallel_1.0.17          
##   [9] lattice_0.21-9              MASS_7.3-60                
##  [11] backports_1.5.0             magrittr_2.0.3             
##  [13] limma_3.58.1                sass_0.4.9                 
##  [15] rmarkdown_2.27              jquerylib_0.1.4            
##  [17] yaml_2.3.9                  plotrix_3.8-4              
##  [19] httpuv_1.6.15               cowplot_1.1.3              
##  [21] RColorBrewer_1.1-3          multcomp_1.4-25            
##  [23] abind_1.4-5                 zlibbioc_1.48.2            
##  [25] GenomicRanges_1.54.1        BiocGenerics_0.48.1        
##  [27] RCurl_1.98-1.14             TH.data_1.1-2              
##  [29] sandwich_3.1-0              circlize_0.4.16            
##  [31] GenomeInfoDbData_1.2.11     IRanges_2.36.0             
##  [33] KMsurv_0.1-5                S4Vectors_0.40.2           
##  [35] codetools_0.2-19            DelayedArray_0.28.0        
##  [37] DT_0.33                     tidyselect_1.2.1           
##  [39] shape_1.4.6.1               farver_2.1.2               
##  [41] matrixStats_1.3.0           stats4_4.3.2               
##  [43] jsonlite_1.8.8              GetoptLong_1.0.5           
##  [45] dr4pl_2.0.0                 survival_3.5-7             
##  [47] iterators_1.0.14            systemfonts_1.1.0          
##  [49] foreach_1.5.2               tools_4.3.2                
##  [51] ragg_1.3.2                  Rcpp_1.0.12                
##  [53] glue_1.7.0                  SparseArray_1.2.4          
##  [55] xfun_0.45                   mgcv_1.9-0                 
##  [57] MatrixGenerics_1.14.0       piano_2.18.0               
##  [59] GenomeInfoDb_1.38.8         shinydashboard_0.7.2       
##  [61] withr_3.0.0                 fastmap_1.2.0              
##  [63] fansi_1.0.6                 shinyjs_2.1.0              
##  [65] relations_0.6-13            caTools_1.18.2             
##  [67] digest_0.6.36               timechange_0.3.0           
##  [69] R6_2.5.1                    mime_0.12                  
##  [71] textshaping_0.4.0           colorspace_2.1-0           
##  [73] gtools_3.9.5                tensor_1.5                 
##  [75] utf8_1.2.4                  generics_0.1.3             
##  [77] renv_1.0.7                  data.table_1.15.4          
##  [79] htmlwidgets_1.6.4           S4Arrays_1.2.1             
##  [81] sets_1.0-25                 pkgconfig_2.0.3            
##  [83] gtable_0.3.5                XVector_0.42.0             
##  [85] survMisc_0.5.6              htmltools_0.5.8.1          
##  [87] carData_3.0-5               fgsea_1.28.0               
##  [89] clue_0.3-65                 scales_1.3.0               
##  [91] Biobase_2.62.0              png_0.1-8                  
##  [93] km.ci_0.5-6                 rstudioapi_0.16.0          
##  [95] tzdb_0.4.0                  rjson_0.2.21               
##  [97] nlme_3.1-163                visNetwork_2.1.2           
##  [99] cachem_1.1.0                zoo_1.8-12                 
## [101] GlobalOptions_0.1.2         KernSmooth_2.23-22         
## [103] pillar_1.9.0                vctrs_0.6.5                
## [105] maxstat_0.7-25              gplots_3.1.3.1             
## [107] slam_0.1-50                 promises_1.3.0             
## [109] ggpubr_0.6.0                car_3.1-2                  
## [111] xtable_1.8-4                cluster_2.1.4              
## [113] evaluate_0.24.0             mvtnorm_1.2-5              
## [115] cli_3.6.3                   compiler_4.3.2             
## [117] rlang_1.1.4                 crayon_1.5.3               
## [119] ggsignif_0.6.4              survminer_0.4.9            
## [121] labeling_0.4.3              marray_1.80.0              
## [123] stringi_1.8.4               BiocParallel_1.36.0        
## [125] munsell_0.5.1               Matrix_1.6-1.1             
## [127] hms_1.1.3                   statmod_1.5.0              
## [129] shiny_1.8.1.1               SummarizedExperiment_1.32.0
## [131] highr_0.11                  rbibutils_2.2.16           
## [133] drc_3.0-1                   exactRankTests_0.8-35      
## [135] igraph_2.0.3                broom_1.0.6                
## [137] bslib_0.7.0                 fastmatch_1.1-4